########################################################################
# PLS 206 — Week 5 live coding, part 2: explore, then check
# Grey Monroe, Fall 2026
########################################################################

# install.packages(c("ggplot2", "MuMIn", "car", "MASS"))
library(ggplot2)
library(MuMIn)
library(car)
library(MASS)

########################################################################
# 1. Pistacia: bio7, VIF
########################################################################

pis <- read.csv("https://greymonroe.github.io/PLS_206/data/pistacia_niche.csv")
names(pis) <- sub("^(bio[0-9]+)_.*", "\\1", names(pis))   # bio1 ... bio19
pis$y <- as.integer(pis$species_short == "terebinthus")    # 1 = P. terebinthus
table(pis$species_short)

# a first model: species ~ annual mean temperature
m_1 <- glm(y ~ bio1, family = binomial, data = pis)
p_ter <- fitted(m_1)                                   # P(terebinthus), each record
sum(log(ifelse(pis$y == 1, p_ter, 1 - p_ter)))         # likelihood by hand
logLik(m_1)                                            # the same

# all 19 at once
bios <- paste0("bio", 1:19)
m_all <- glm(y ~ ., family = binomial, data = pis[, c("y", bios)])
round(coef(m_all)[c("bio5", "bio6", "bio7")])   # +-30,000: nonsense
range(pis$bio7 - (pis$bio5 - pis$bio6))         # bio7 = bio5 - bio6

# sampling effort (log10), complete rows, climate in SD units
pis$effort <- log10(pis$plant_obs_effort)
pis <- subset(pis, !is.na(effort))
pis[bios] <- scale(pis[bios])

# drop the worst VIF, refit, repeat
keep <- setdiff(bios, "bio7")
repeat {
  v <- vif(glm(reformulate(keep, "y"), family = binomial, data = pis))
  if (max(v) < 5) break
  keep <- setdiff(keep, names(which.max(v)))
}
keep                            # 7 left
round(v, 1)

# what VIF means: predict one variable from the others
r2 <- summary(lm(bio6 ~ bio8 + bio4 + bio10 + bio15 + bio9 + bio2 + bio12 + bio3, data = pis))$r.squared
r2                              # 0.989: bio6 is almost a mix of the other 8
1 / (1 - r2)                    # its VIF

########################################################################
# 2. dredge(), importance, averaging
########################################################################

options(na.action = "na.fail")  # dredge() insists
full <- glm(y ~ bio1 + bio3 + bio4 + bio8 + bio9 + bio18 + bio19 + effort,
            family = binomial, data = pis)
dd <- dredge(full, fixed = ~ effort)   # effort in every model
nrow(dd)                        # 128
dd[1:10, recalc.weights = FALSE]       # top model: 37% of the weight
subset(dd, delta < 2, recalc.weights = FALSE)

dredge(full, fixed = ~ effort, rank = "BIC")[1:5, recalc.weights = FALSE]   # smaller top model

sw(dd)                          # importance: summed weights

avg <- model.avg(get.models(dd, subset = delta < 2))
round(coef(avg, full = TRUE), 2)

# the averaged model against the data: bio4, others at their mean (0)
nd <- data.frame(matrix(0, 100, length(keep), dimnames = list(NULL, keep)))
nd$effort <- mean(pis$effort)
nd$bio4 <- seq(-2.5, 2.5, length.out = 100)
pr <- predict(avg, nd, se.fit = TRUE, full = TRUE)
nd$fit <- plogis(pr$fit)
nd$lwr <- plogis(pr$fit - 1.96 * pr$se.fit)
nd$upr <- plogis(pr$fit + 1.96 * pr$se.fit)

p <- ggplot(pis, aes(bio4, y)) +
  geom_jitter(height = 0.03, width = 0, alpha = 0.05) +
  labs(x = "Temperature seasonality, bio4 (SD units)", y = "P(P. terebinthus)") +
  theme_classic()
p
p + geom_ribbon(data = nd, aes(bio4, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = nd, aes(bio4, fit), colour = "blue")

########################################################################
# 3. stepAIC on resampled farms
########################################################################

farms <- read.csv("https://greymonroe.github.io/PLS_206/data/cover_crop.csv")
farms <- subset(farms, select = -c(farm, clay))
P <- setdiff(names(farms), "biomass")

step_f <- stepAIC(lm(biomass ~ ., data = farms), trace = 0)
step_f$anova                          # the path: one predictor dropped per step
attr(terms(step_f), "term.labels")    # gdd + tmean_winter; the truth is gdd alone

# forward: start with none, add the best one each step
step_fw <- stepAIC(lm(biomass ~ 1, data = farms), direction = "forward", trace = 0,
                   scope = list(lower = ~ 1, upper = reformulate(P)))
step_fw$anova                         # + gdd, + tmean_winter, stop

# collect the data again, 200 times
set.seed(2)
sel <- replicate(200, {
  s <- farms[sample(50, 50, replace = TRUE), ]
  m <- stepAIC(lm(biomass ~ ., data = s), trace = 0)
  P %in% attr(terms(m), "term.labels")
})
rownames(sel) <- P
sort(rowMeans(sel), decreasing = TRUE)   # gdd 99%; junk 35-58%
table(colSums(sel))                      # how many it kept: 2 to 13

########################################################################
# 4. Pure noise
########################################################################

set.seed(41)
noise <- data.frame(y = rnorm(100), matrix(rnorm(100 * 20), 100))
step_n <- stepAIC(lm(y ~ ., data = noise), trace = 0)
round(summary(step_n)$coefficients, 3)   # p = 0.006 and 0.001, from random numbers

# 1,000 times: choose on 50 rows, test on the same 50 or on the other 50
set.seed(1)
fz <- replicate(1000, {
  nz <- data.frame(y = rnorm(100), matrix(rnorm(100 * 20), 100))
  a <- nz[1:50, ]; b <- nz[51:100, ]
  m <- stepAIC(lm(y ~ ., data = a), trace = 0)
  tl <- attr(terms(m), "term.labels")
  if (length(tl) == 0) return(c(same = FALSE, other = FALSE))
  c(same  = any(summary(m)$coefficients[-1, 4] < 0.05),
    other = any(summary(lm(reformulate(tl, "y"), data = b))$coefficients[-1, 4] < 0.05))
})
rowMeans(fz)                             # 86% vs 25%

########################################################################
# 5. AIC and leave-one-out prediction
########################################################################

# add predictors best first; for each model, AIC and leave-one-out error
ord <- names(sort(-abs(cor(farms)[P, "biomass"])))
res <- t(sapply(0:15, function(k) {
  m <- lm(if (k == 0) biomass ~ 1 else reformulate(ord[1:k], "biomass"), data = farms)
  loo <- sqrt(mean((resid(m) / (1 - hatvalues(m)))^2))   # predict each farm from the other 49
  c(k = k, AIC = AIC(m), LOO = loo)
}))
res
res[which.min(res[, "AIC"]), "k"]    # 2
res[which.min(res[, "LOO"]), "k"]    # 2: the same